% script to process IPM and PPT data

% The data are 1 month 1 sec for Feb 2010
ipm = load('/Users/manojnair/projects/tsunami/IPM');
ppt = load('/Users/manojnair/projects/tsunami/PPT');
% selected filtered data
ipm = load('/Users/manojnair/projects/tsunami/solartides/output.txt');

L = ppt > 99999.0;
ppt(L) = NaN;


L = ipm > 99999.0;
ipm(L) = NaN;

L = isnan(ipm(:,1)) | isnan(ipm(:,2)) | isnan(ipm(:,3));

xaxis = 1:length(ipm);


ipm(:,1) = interp1(xaxis(~L), ipm(~L,1), xaxis, 'linear');
ipm(:,2) = interp1(xaxis(~L), ipm(~L,2), xaxis, 'linear');
ipm(:,3) = interp1(xaxis(~L), ipm(~L,3), xaxis, 'linear');


L = isnan(ppt(:,1)) | isnan(ppt(:,2)) | isnan(ppt(:,3));

xaxis = 1:length(ppt);


ppt(:,1) = interp1(xaxis(~L), ppt(~L,1), xaxis, 'linear');
ppt(:,2) = interp1(xaxis(~L), ppt(~L,2), xaxis, 'linear');
ppt(:,3) = interp1(xaxis(~L), ppt(~L,3), xaxis, 'linear');


% save('/Users/manojnair/projects/tsunami/ipm.txt', 'ipm','-ascii',  '-tabs' );
% save('/Users/manojnair/projects/tsunami/ppt.txt', 'ppt','-ascii',  '-tabs' );

% plot the data

fday = linspace(datenum(2010,02,1,0,0,0),datenum(2010,02,28,23,59,59),28*24*3600);
subplot(311);
plot(fday,ipm(:,1)-nanmean(ipm(:,1)),'b');
set(gca,'FontSize',16);
datetick('x','keeplimits');
hold on
plot(fday,ppt(:,1)-nanmean(ppt(:,1)),'r');
ylabel('\Delta X nT');
subplot(312);
plot(fday,ipm(:,2)-nanmean(ipm(:,2)),'b');
set(gca,'FontSize',16);
datetick('x','keeplimits');
hold on
plot(fday,ppt(:,2)-nanmean(ppt(:,2)),'r');
ylabel('\Delta Y nT');

subplot(313);
plot(fday,ipm(:,3)-nanmean(ipm(:,3)),'b');
set(gca,'FontSize',16);
datetick('x','keeplimits');
hold on
plot(fday,ppt(:,3)-nanmean(ppt(:,3)),'r');
ylabel('\Delta Z nT');


%% Univariate TF

[Txy,   F] = tfestimate(ppt(:,1),ipm(:,1),hanning(2048*10),1,2048*10,1);
[Txy_1, F] = tfestimate(ppt(:,2),ipm(:,2),hanning(2048),1,2048,1);
[Txy_2, F] = tfestimate(ppt(:,1),ipm(:,3),hanning(2048),1,2048,1);
[Txy_3, F] = tfestimate(ppt(:,3),ipm(:,1),hanning(2048),1,2048,1);


[Cxy,   F] = mscohere(ppt(:,1),ipm(:,1),hanning(2048),1,2048,1);
[Cxy_1, F] = mscohere(ppt(:,2),ipm(:,2),hanning(2048),1,2048,1);
[Cxy_2, F] = mscohere(ppt(:,1),ipm(:,3),hanning(2048),1,2048,1);
[Cxy_3, F] = mscohere(ppt(:,3),ipm(:,1),hanning(2048),1,2048,1);


% plot the data

loglog((1./F), abs(Txy),'r-', 'LineWidth', 2);
hold on;
set(gca, 'FontSize', 16);
loglog((1./F), abs(Txy_1),'b-', 'LineWidth', 2);
loglog((1./F), abs(Txy_2),'k-', 'LineWidth', 2);
loglog((1./F), abs(Txy_3),'g-', 'LineWidth', 2);
legend('PPTX-IPMX','PPTY-IPMY','PPTX-IPMZ','PPTZ-IPMX');
xlabel('Period in S');
ylabel('TF Magnitude');

semilogx((1./F), angle(Txy) * 180/pi,'r-', 'LineWidth', 2);
hold on;
set(gca, 'FontSize', 16);
semilogx((1./F), angle(Txy_1) * 180/pi,'b-', 'LineWidth', 2);
semilogx((1./F), unwrap(angle(Txy_2)) * 180/pi,'k-', 'LineWidth', 2);
semilogx((1./F), unwrap(angle(Txy_3)) * 180/pi,'g-', 'LineWidth', 2);
legend('PPTX-IPMX','PPTY-IPMY','PPTX-IPMZ','PPTZ-IPMX');
xlabel('Period in S');
ylabel('TF Phase');


%% Bivariate TF


nfft = 4096 * 5 ;

[HxHx,   F] = pwelch(ppt(:,1),hanning(nfft),1,nfft,1);
[HyHy,   F] = pwelch(ppt(:,2),hanning(nfft),1,nfft,1);
[HxHy,   F] = cpsd(ppt(:,1),ppt(:,2),hanning(nfft),1,nfft,1);
[HzHx,   F] = cpsd(ppt(:,3),ppt(:,1),hanning(nfft),1,nfft,1);
[HzHy,   F] = cpsd(ppt(:,3),ppt(:,2),hanning(nfft),1,nfft,1);
[HyHx,   F] = cpsd(ppt(:,2),ppt(:,1),hanning(nfft),1,nfft,1);

% Txy = (HxHx .* (HzHy) - HxHy .* (HzHx))./(HxHx .* (HyHy) - HxHy .* (HyHx));
% 
% subplot(211);
% loglog((1./F), abs(Txy),'b-', 'LineWidth', 2);
% axis([ 50       10000       -inf 2]);
% 
% subplot(212);
% semilogx((1./F), unwrap(angle(Txy)) * 180/pi,'b-', 'LineWidth', 2);
% axis([ 50       10000       -180 180]);
% 
% 

for i = 1:length(HzHx);
    
EH=[HzHx(i) HzHy(i)];

HH=[HxHx(i) HxHy(i); HyHx(i) HyHy(i)];

Z(i,:) = EH/HH;

end;

%% Plot

L = isinf((1./F)) | (1./F) < 50;

subplot(211);
loglog((1./F(~L)), abs(Z(~L,2)),'b*', 'LineWidth', 1);
axis([ 50       10000       -inf 2]);
grid on
hold on
loglog((1./F(~L)), abs(Z(~L,1)),'r*', 'LineWidth', 1);
subplot(212);
semilogx((1./F(~L)), angle(Z(~L,2))* 180/pi,'b*', 'LineWidth', 1);
hold on
semilogx((1./F(~L)), angle(Z(~L,1))* 180/pi,'r*', 'LineWidth', 1);
axis([ 50       10000       -180 180]);

grid on


subplot(211);
set(gca, 'FontSize', 16);
%legend('TF1 (Matlab)','TF2 (Matlab)','TF1 (EMTF)','TF2 (EMTF)');
legend('TF1','TF2');

xlabel('Period (S)');
ylabel(' TF Magnitude');
subplot(212);
set(gca, 'FontSize', 16);
%legend('TF1 (Matlab)','TF2 (Matlab)','TF1 (EMTF)','TF2 (EMTF)');
legend('TF1','TF2');

xlabel('Period (S)');
ylabel(' TF Phase');
subplot(211);
%title('Bx_{ppt} = TF1*Bx_{ipm} + TF2*By_{ipm}');
%title('Bz_{IPM} = TF1*Bx_{IPM} + TF2*By_{IPM}');
%title('Bx_{IPM} = TF1*Bx_{PPT} + TF2*By_{PPT}')
%title('Bz_{PPT} = TF1*Bx_{IPM} + TF2*By_{IPM}')
%title('Bz_{IPM} = TF1*Bx_{PPT} + TF2*By_{PPT}')
title('Bz_{PPT} = TF1*Bx_{PPT} + TF2*By_{PPT}')


subplot(212);
%title('Bx_{ppt} = TF1*Bx_{ipm} + TF2*By_{ipm}');
%title('Bz_{IPM} = TF1*Bx_{IPM} + TF2*By_{IPM}');
%title('Bz_{IPM} = TF1*Bx_{PPT} + TF2*By_{PPT}')
%title('Bz_{PPT} = TF1*Bx_{IPM} + TF2*By_{IPM}')
%title('Bz_{IPM} = TF1*Bx_{PPT} + TF2*By_{PPT}')
title('Bz_{PPT} = TF1*Bx_{PPT} + TF2*By_{PPT}')





%% Get the TF from Egbert's codes

%cfile= '/Users/manojnair/projects/tsunami/x1y1z1z1z2.zss'
cfile= '/Users/manojnair/projects/tsunami/x2y2z2z2z1.zss'

TF = read_egbert_tf(cfile);
% 1=ipm and 2= ppt

% plot the data

Z_emtf = TF.TF;
periods = TF.T;
err = TF.SIG_s;

subplot(211);

loglog(periods, abs(Z_emtf(3,:)),'r.-');

hold on;
 for i = 1: length(Z_emtf);
    plot([periods(i) periods(i)], [abs(Z_emtf(1,i)) - ... 
        abs(err(1,i))/abs(Z_emtf(1,i)) abs(Z_emtf(1,i)) + abs(err(1,i))/abs(Z_emtf(1,i))],'k','LineWidth',1);
    hold on;
end;

hold on;
loglog(periods, abs(Z_emtf(4,1:25)),'b.-');
% loglog(periods, abs(Z(3,1:25)),'k');
% loglog(periods, abs(Z(4,1:25)),'c');
% legend('TF1','TF2','TF3','TF4');

subplot(212);

semilogx(periods, angle(Z_emtf(3,1:25)) * 180/pi,'r.-');
hold on;
semilogx(periods, angle(Z_emtf(4,1:25)) * 180/pi,'b.-');
% semilogx(periods, angle(Z(3,1:25)) * 180/pi,'k');
% semilogx(periods, angle(Z(4,1:25)) * 180/pi,'c');
% legend('TF1','TF2','TF3','TF4');
